Time-Dependent Hartree-Fock Approach to Nuclear Pasta at Finite Temperature 
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We present simulations of neutron-rich matter at subnuclear densities, like supernova matter, with 
the time-dependent Hartree-Fock approximation at temperatures of several MeV. The initial state 
consists of a particles randomly distributed in space that have a Maxwell-Boltzmann distribution 
in momentum space. Adding a neutron background initialized with Fermi distributed plane waves 
the calculations reflect a reasonable approximation of astrophysical matter. This matter evolves 
into spherical, rod-like, and slab-like shapes and mixtures thereof. The simulations employ a full 
Skyrme interaction in a periodic three-dimensional grid. By an improved morphological analysis 
based on Minkowski functionals, all eight pasta shapes can be uniquely identified by the sign of only 
two valuations, namely the Euler characteristic and the integral mean curvature. 

PACS numbers: 26.50.+x,21.60. Jz,21.65.-f,97.60.BW 



I. INTRODUCTION 

Core-collapse supernovae (see, 0, Q for a review) are 
considered to be relevant for many astrophysical phenom- 
ena including synthesis of heavy nuclei, acceleration of 
cosmic rays, and formation of neutron stars and black 
holes. The mechanism of such supernovae, however, has 
yet to be clarified, although hydrodynamics simulation 
studies of stellar collapse and subsequent explosion con- 
tinue to be improved with respect to dimensionality and 
neutrino transport as well as inclusion of general rela- 
tivity, rotation, and magnetic fields Q. For such clar- 
ification, the properties of matter that constitutes the 
supernova cores such as the equation of state, the nu- 
clei present, and the neutrino opacity are indispensable. 
This matter, which is believed to be charge-neutral and 
mostly beta-equilibrated while having a degenerate gas 
of electron neutrinos trapped inside, is often referred to 
as supernova matter. In early studies of supernova mat- 
ter below normal nuclear density and at temperatures of 
order or lower than 10 MeV, nuclei in equilibrium were 
shown to be relatively neutron-rich by using a liquid-drop 
model Q , a Thomas-Fermi theory [j| , and a Hartree-Fock 
(HF) theory @. 

At densities just below normal nuclear density, nearly 
spherical nuclei arranged in a Coulomb lattice and em- 
bedded in a roughly uniform neutralizing background of 
electrons are so closely packed that the total surface area 
becomes very large. Then, the system tends to lower 
the total surface area by elongating the nuclei and form- 
ing nuclear rods 0, Hj]. It was predicted from seminal 
liquid-drop calculations 0, [l(| that with further increase 
in density, the shapes of inhomogeneous nuclear matter 
change from rods to slabs, tubes, and bubbles until the 
system melts into uniform matter. These exotic nuclear 
configurations are often denoted as pasta nuclei, which 
arise from subtle competition between the surface and 
Coulomb energies. The pasta nuclei are manifestations of 
liquid-gas mixed phases of nuclear matter, among which 



the phases with rods, slabs, and tubes can be regarded 
as liquid crystals flTj . 

In the conventional liquid-drop approach 0, [l2| , pasta 
nuclei were studied by assuming a geometrical shape 
of nuclear matter and using the Wigner-Seitz approxi- 
mation. For more realistic description of pasta nuclei, 
multi-dimensional calculations without any assumption 
of the geometrical shape were performed in the Thomas- 
Fermi approach [lU, [lj] , the Hartree-Fock approach [lj| , 
and the quantum molecular dynamics (QMD) approach 
[l6l lvi\ . In these calculations, the shapes assumed in 
the Wigner-Seitz approximation were reproduced, while 
some cases did indicate the presence of more complex 
structures. 

Astrophysically, the possible influence of pasta nuclei 
on the neutrino opacity may affect supernova explosions 
and subsequent protoncutron star cooling. The cross sec- 
tion for the neutrino-nucleus scattering depends on the 
structure of inhomogeneous nuclear matter fl8l - l20j . In 
fact, the neutrino scattering processes are no longer co- 
herent in the directions in which nonspherical nuclei are 
elongated. This is a great contrast to the case of roughly 
spherical nuclei of which the finiteness in any direction 
leads to constructive interference in the scattering. 

Recently, a time-dependent Hartree-Fock (TDHF) ap- 
proach has been utilized to describe pasta nuclei in zero- 
temperature supernova matter pH [22j . In contrast to 
static calculations which concentrate on the one config- 
uration at minimal energy, the TDHF approach allows 
nucleons to explore multiple low-lying configurations in 
the energy landscape mapped by the time evolution of 
the mean field. In this work we perform such dynamic 
calculations at finite temperature, which are expected to 
give a realistic description of supernova matter as long 
as two-nucleon collisions can be ignored. As we shall see, 
the resultant map of pasta shapes is basically consistent 
with the phase diagrams obtained from the QMD calcu- 
lations jl7| . 

The TDHF simulations produce a great variety of very 
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involved, three dimensional structures. One needs tools 
to characterize a given configuration in terms of a few 
key numbers. This is simple for single finite nuclei where 
radius and a couple of deformation parameters serve very 
well. However, we encounter here compounds where 
much less is known ahead of time. For example, we do 
not even know the number of subunits or their connectiv- 
ity. A very general and powerful means to quantify com- 
plex structures is provided by the Minkowski functionals, 
shape descriptors from integral geometry [23|, [24| . These 
measures of structure are important tools in statistical 
physics (25l - [27| and pattern analysis [H, [29| . They were 
already introduced to astronomy for structure analysis of 
galaxy clusters in Ref. and later used to investigate 
both point processes in cosmology [HI H2| and the cosmic 
microwave background [331 ] . Such Minkowski function- 



als provide, of course, the optimal technique to quantify 
pasta matter and so classify the different possible shapes. 
This had already been exploited in [34H36I ] . We aim to 
continue and deepen studies along this line. 

In Sec. II, the TDHF approach is briefly reviewed. In 
Sec. Ill, a finite-temperature map of pasta shapes for in- 
homogeneous nuclear matter is presented. With the help 
of Minkowski functionals given in Sec. IV, we do mor- 
phological identifications of various phases in the map of 
pasta shapes in Sec. V. In Sec. VI we investigate thermal 
fluctuations in the density profiles and discuss their con- 
nection with melting of the inhomogeneous structures. 
Concluding remarks are given in Sec. VII. 



II. TDHF DESCRIPTION 

In this framework, we want to describe the pasta struc- 
ture with the time-dependent Hartree-Fock (TDHF) ap- 
proximation. TDHF is a self-consistent mean field the- 
ory which was originally proposed by Dirac in 1930 [37j . 
In general, the TDHF equations are derived from a time- 
dependent variational principle with the variational space 
restricted to one single, time dependent Slater determi- 
nant. In nuclear applications, TDHF usually employs an 
effective interaction as, e.g., the Skyrme force, for a re- 
view see [38j . The treatment becomes then practically 
identical to a time-dependent density functional theory 
in local density approximation as it is widely used in elec- 
tronic systems, see e.g. [39]. Computational limitations 
restricted earlier TDHF calculations in system size or 
symmetry |4l[ . With the advance of computational 
power, it is now possible to perform full 3D calculations 
without any symmetry restrictions even for a large num- 
ber of wave functions. This led to a revival of nuclear 
TDHF studies in various dynamical regimes as, e.g., col- 
lective vibrations [HI, [43| , large amplitude motion in fu- 
sion and nuclear reactions Q [4f| , or, similar as done 
here, for simulations of stellar matter [2ll [22j. 

The TDHF calculations are performed on an equidis- 
tant 3D grid in coordinate space with 16 grid points in 
each direction with a grid spacing of 1 fm. As usually 



done in the physics of extended systems, we simulate 
infinite matter in a finite box with periodic boundary 
conditions. The Coulomb problem is solved with as- 
sumed global neutrality which means that we compen- 
sate the proton charges by a homogeneous negatively 
charged background which is supposed to simulate the 
electron gas around. This amounts to use the Coulomb 
solver with periodic boundary conditions and skipping 
the zero momentum component of the Coulomb field in 
Fourier space. The kinetic energy operator is evaluated 
in Fourier space and time evolution is performed using a 
Taylor expansion of the mean-field propagator. For the 
present calculations we choose Skyrme parametrization 
SLy6 which was developed with an emphasis on describ- 
ing neutron rich matter (46j. 



III. MAP OF REGIMES 
A. Initialization 

We consider matter with a proton fraction of X p = 1 /3 . 
In order to produce not too high excitations, we initialize 
the calculations by distributing a number N a of a parti- 
cles randomly over the grid keeping a minimal distance 
between each a core to avoid perturbations by the over- 
lap. To that end, we use ground wavef unctions for the 
a particles produced before in a stationary HF calcula- 
tion. The number N a is varied to produce systems with 
different densities. These a cores are distributed in mo- 
mentum space following a Maxwell-Boltzmann distribu- 
tion with a given initialization temperature T- m i t . Further 
neutron states are introduced as plane wave states using 
a Fermi distribution with the same temperature to de- 
scribe free background neutrons. Note that always both 
spin up and spin down states are occupied simultaneously 
to avoid a large spin excitation. Finally, before starting 
the TDHF calculation, these wave functions have to be 
ortho- normalized. 

This randomized setup leads usually to a system with- 
out any symmetry. The time evolution is assumed to be 
sufficiently chaotic such that robust final states proper- 
ties are reached which depend only on temperature and 
average density, although there are some regions where 
two states compete and the final outcome depends also 
on the initialization. The average density is tuned by 
the number of particles and the temperature by Tinit in 
the sampling of a particles and neutrons. Even with 
Tinit = we store a considerable amount of excitation in 
the system stemming from the interaction energy of the 
a particles. 



B. Temperature 

The excitation state of matter is characterized by its 
actual temperature T. It is, however, difficult to mea- 
sure temperature in the quantum simulation of a micro- 
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FIG. 1. Typical shapes of pasta structures at the lowest tem- 
peratures. These continue to evolve with time but do not 
change their morphological character anymore. Bubble shape 
illustrations show gas phase, which is indicated by the color- 
scale. 



canonical ensemble as generated by the present TDHF 
calculations. As a rough measure, we estimate T by 
the excitation energy E* through the relation E* = 
(tt 2 A)(2ef) T 2 where A is the total number of nuclcons 
and £p the Fermi energy of the system [47J . The excita- 
tion energy E* is defined as the difference of the actual 
energy E to the ground state energy Eq. The latter is 
computed by solving the static HF problem for the given 
density. This amounts to 



T 



2s F 



E* 



(1) 



This excitation temperature T is to be distinguished 
from the "initial temperature" Tj n i t which we use to 
boost stochastically the initial ensemble of a particles 
and the background neutrons. Calculating the excita- 
tion temperature T for the different values of the ini- 
ti al temperature Tjnit, w e find empirically T « 7 MeV + 
\/Tinit(Ti„it + 100 MeV)/6. This is taken henceforth for 
a rough calibration of T. 




0.06 0.08 0.1 0.12 

P [fm- 3 ] 

□ sphere Drod Drod(2) ■rod(3) ■slab 

□ rod(2) bubble Brod bubble dsphere bubble □ uniform 

FIG. 2. Map of pasta shapes achieved in a TDHF calculation 
starting from a gas of a particles with neutron background for 
various real temperatures and mean densities. The proton 
fraction is 1/3. Each dot represents two calculations. The 
solid black line shows a phase separation line discussed in 
Sec. EH 



strong fluctuations can be observed. 

In FIG. [1] many different pasta shapes arc classified. 
Among the structures found are rod and slab struc- 
tures which have been discovered, e.g., by QMD calcu- 
lations (3(1 . Rod(2) corresponds to rods forming a two- 
dimensional layer. The shape rod(3) describes three rods 
in x-, y- and z-direction which cross in one point. Similar 
shapes were discovered in [2l|, |48| and the rod(3) struc- 
ture was also found in (l5j . 

FIG. [2] shows for which density and temperature the 
different shapes appear. Note that two calculations for 
each point in the map were performed. As we have ran- 
dom initial conditions, the final shapes may differ al- 
though the same values for temperature and mean den- 
sity are assumed. The hatched areas with a mixture of 
colors indicate the different final states reached in these 
cases. At higher densities, bubble structures can be seen 
(the gas phase has a pasta like shape). Many shapes can 
coexist here in one point of the map. 

In the region of slab and rod(3) structures, the final 
state which is formed seems to highly depend on the ini- 
tial condition. The similarity of these two structures is, 
that the liquid phase and gas phase have the same shapes. 
So gas and liquid phase symmetrically complement each 
other and fill out almost the same volume. 



IV. MINKOWSKI FUNCTIONALS 



C. First overview 

Each setup is evolved in time for 1500 fm/c. After that 
time, shapes do not change significantly. The system goes 
over into a pasta state where some type of equilibrium 
is achieved. Depending on the temperature differently 



As already outlined in section [IJ a powerful tool to 
quantify the involved pasta shapes are Minkowski func- 
tionals [26|. There are four Minkowski functionals W v 
defined for a spatial domain K in three dimensions: they 
are proportional to its volume Wq oc V, its surface-area 
Wi (x A, the integral mean curvature W2 oc J gK dA (k\ + 
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TABLE I. Minkowski functionals in three dimensions evalu- 
ated for a domain K, with Ki and k 2 the principle curvatures 
on dK. 



Wo 
Wi 

w 2 
w 3 



oc 
cx 
cx 
cx 



V 

Sbk^ 

Sok dA( K i + K2 )/2 
fax dA Kl ' K2 



volume 
surface 

int. mean curvature 
Euler characteristic 



«2)/2, and the topological Eulcr-Poincare characteristic 
W3 oc x, which is equal to the integrated Gaussian cur- 
vature f dK dAni ■ K2- Here, k\ and K2 are the principle 
curvatures on dK, the bounding surface of K. The Euler 
characteristic is a topological constant. TAB. U summa- 
rizes these definitions of the Minkowski functionals. 

Minkowski functionals are useful shape descriptors. 
They are robust against noise and have short computa- 
tion times, both due to their additivity property. These 
shape indices form a complete basis of all functionals de- 
fined on unions of convex sets which are translational and 
rotational invariant, additive, and at least continuous on 
convex sets [49|. Thus, the Minkowski functionals pro- 
vide a complete scalar morphological characterization of 
the pasta shapes these functionals. The Minkowski func- 
tionals can be generalized to Minkowski tensors for an 
anisotropic structure characterization [27| ■ 

The TDHF results for pasta matter are given as gray- 
scale data with a density pi assigned to each voxel i (= 
grid point). In order to compute the Minkowski function- 
als of a domain K , the density field is turned into binary 
data via thresholding. A threshold density pth is intro- 
duced. Each voxel i with pi > p t h is set to black, but if 
Pi < Pth, white is assigned. The Minkowski functionals 
of the union of all black voxels, interpreted as a polygon 
K, can quantify the shape of the pasta matter in depen- 
dence of Qjh. Beforehand however, the marching cube al- 
gorithm [281 ] is applied, which provides a smoothed polyg- 
onal representation in order to reduce voxelization errors. 
The Minkowski functionals of the polygon are evaluated, 
e.g. , using the linear-time algorithm from Refs. [U IHoj 

m- 

The threshold density can take on values pth £ 
[0, Pmax]- But for small thresholds p t h ~ and for 
high thresholds pth ~ Pmax, the Minkowski functionals 
show erratic behavior; thus, the structure values for these 
thresholds are ignored [13, [H| . 



V. STRUCTURE CLASSIFICATION 

Calculating the values of the Minkowski scalars for all 
shapes as a function of threshold density (FIG. [3|), we 
see that we can uniquely classify the different shapes in 
simple fashion: only the signs of integral mean curvature 
and Euler characteristic are needed (see TAB. HT|) . 

The central regime with solid lines in FIG. [3] indi- 
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FIG. 3. W2, proportional to the integral mean curvature, and 
W3, proportional to the Euler characteristic as function of 
threshold density for the samples shown in FIG[T] 



cate the physical values for the threshold density which 
yield reliably stable values of W2 and W3 for pasta 
configurations with the lowest value of temperature are 
pth £ [0.03 fm~ 3 , 0.08 fm~ 3 ]. For lower threshold den- 
sities only few dots are recognized as "white" (below 
threshold density) due to quantum fluctuations in the 
gas phase. At higher densities quantum fluctuations of 
the liquid phase are observed. 

Rod(2), rod(2) bubble, and rod(3) consist in channels 
and tunnels; thuSjthe domain has a negative Euler char- 
acteristic x < |53l ]. 

The bubble shapes are the complements of the accord- 
ing sphere, rod, or rod(2) shapes. As the sign of both 
main curvatures Ki and K2 changes, the mean curvature 
H = («i + k 2 )/2 also changes sign and so docs Wi- How- 
ever, the Gaussian curvature G = «i • K2 and thus the 
Euler characteristic remain constant. 

Both the slab and the rod(3) shape are symmetric in 
gas and liquid phase, as described in Sec. Mil The shape 
of the liquid phase at low thresholds corresponds to the 
gas phase at high thresholds. Thus, the integral mean 
curvature as a function of the threshold pth is point sym- 
metric w.r.t. a mid value for the threshold density. This 



TABLE II. Signs of the integral mean curvature and Euler 
characteristic for each observed shape. 
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Pth [ fm " 3 l Pth [ fm " 3 l Pth [ fm " 3 l 

FIG. 4. Threshold density profiles for calculations varying 
in density and temperature. On each single plot the number 
of grid points n corresponding to a certain density value is 
plotted. The dashed lines show the mean density. 



is trivially true for the slab with W2 ~ 0, but holds also 
for the non-zero mean curvature of rod(3). 

The results indicate that for a complete classification 
of the observed pasta shapes, the integral mean curvature 
function of pth must be known. 

Especially at higher temperatures and mean densities 
mixtures of phases can appear, in a sense that at different 
threshold densities the pasta matter can take on different 
shapes (cf. [HI])- With increasing threshold densities 
the shapes can change to shapes which appear usually at 
lower mean densities e.g. a rod(2) bubble shape can take 
on a slab shape at higher values of pth- For these cases 
we take the pasta phase with the bigger range in p t h ■ 



VI. FLUCTUATIONS 

To show the fraction of liquid and gas phase and the 
transition to uniform matter, we plotted the threshold 
density profiles for various calculations in FIG. |U For 
the lowest temperature and low mean density, there is a 
big peak at small densities and a long tail to high densi- 
ties. With increasing mean density a peak at high density 
develops. At the densities where rod(3) and slab shapes 
appear (5 ■ lCT 2 fm -3 < p < 7 ■ 10~ 2 fm~ 3 ) the peaks 
for low and high densities have nearly the same width 
and height. At higher mean densities the peak for low 
densities disappears and only a tail remains. 

With increasing temperature the double peak struc- 
ture vanishes and finally a single peak around the mean 



density is forming which moves successively to central 
density with ever higher T. Consider, e.g., the case of 
the mean density of 6.4 • 10~ 2 fm~ 3 in FIG. [U We can 
observe a double peak structure for the lowest temper- 
ature T = 7MeV. For T = 11.2 Me V we can separate 
the two phases at a mid value of the peak and still ob- 
serve pasta structure. For higher temperatures no pasta 
structure can be observed. Only fluctuations around the 
mean value plotted as a dashed line remain, so that it is 
considered as uniform matter. 

As the pasta structures vanish with higher density, the 
variance from the mean value of the calculation of these 
plots decreases. Going through the systematics of the 
results, we found as a reasonable value for a limit of 
the variance to observe pasta structures 10 -4 fm -6 . The 
phase separation line computed with this observable is 
displayed in FIG. [2] as a bold line. 



VII. CONCLUSION 

In this work we have studied the nuclear composition 
of supernova matter, i.e. matter at sub-nuclear densi- 
ties. To this end, we used self-consistent nuclear mean- 
field theory at the level of time-dependent Hartree-Fock 
(TDHF) using the Skyrme energy functional. We solve 
the TDHF equations without any symmetry restrictions 
on a fully three dimensional coordinate-space grid using 
periodic boundary conditions to simulate infinite matter. 
The initial state is prepared by placing N a a particles 
stochastically over the simulation box and adding 2N a 
neutrons in plane waves which amounts to a 2:1 mix of 
neutrons to protons. Scanning cases with different N a 
allows to tunc different average densities. These initial 
states carry already a minimum of excitation energy be- 
cause of the interaction energy of the a particles and 
the neutrons. We have additionally varied the excita- 
tion state by initializing the a particles and the neutrons 
with a certain amount of kinetic energy. As a measure 
for the actual excitation, we introduce a rough estimate 
for the resulting temperature based on the excitation 
energy (relative to the true ground state of the given 
case). The thus given initial states evolve to a topologi- 
cally stable state of pasta matter within about 1000 fm/c. 
This allows us to deduce a map of pasta shapes in the 
plane of temperature and density. Those TDHF results 
agree qualitatively with the phase diagram calculated in 
QMD [36fl . In realm of "rod" structure, we include fur- 
ther sub-structures, coined rod(2) and rod(3), which ap- 
peared in earlier Hartree-Fock calculations [HI, [2l| . 

As a further means for characterizing these rather in- 
volved three-dimensional structures, we use Minkoswki 
functionals which have been proven to be useful mea- 
sures for structures of fuzzy material. As we deal with 
smooth density distributions which cover continuously 
all values between zero and maximum density, we pro- 
duce to 1 step functions by setting a threshold density 
and wc consider the Minkowski functionals as function 
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of this threshold density. This method of analysis allows 
to characterize the various shapes in terms of a few key 
numbers. This allows, e.g., to discriminate reliably the 
different rod structures, or to work out the smooth tran- 
sition to uniform matter at high temperature. By taking 
the variance of the density profile as a complementing 
further observable, we can define a phase border between 
uniform matter and pasta shapes. 
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